/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

// forceNbdMain.c -- main loop involved in the calculation of nonbonded ("Nbd") interactions. 

//#define DEBUG // remove comment to turn on debugging functionalities

#include "dataStructs.h"

extern const VecR *siteCoords, *siteOrientVec;
extern VecR *siteForce, *siteTorque;
extern const RMat *siteRotationMatrix;
extern const Site *site;  
extern const Atom *atom;
extern const VecR region;
extern const real rCutSoluteElse, rCutMax, rCutSolute, rCutWatWat, rCutLipLip;
extern const int *interactionType, *nebrTab, nebrTabLen, nAtoms;
extern const int stepCount, stepLpp, nSites, lpp;
extern const real **eps, **sigSquared, **productElectroMoments;
extern const real sigWat, epsWat, watDipole, hBondWatWat;

/* external function declarations: */
void CalcChgChgInteraction_SF( real*, VecR*, real, real, real, VecR );
void CalcDipDipInteraction( real*, VecR*, VecR*, VecR*, real, VecR, VecR, VecR, real );
void CalcChgDipInteraction_SF( real*, VecR*, VecR*, real, VecR, real, VecR, real );
void CalcDipDipInteraction_SF( real*, VecR*, VecR*, VecR*, real, VecR, VecR, real, VecR, real );
void CalcDipDipInteraction_CS( real*, VecR*, VecR*, VecR*, real, VecR, VecR, real, VecR, real );
void CalcLennJonesInteraction( real*, VecR*, real, real, real, VecR );
void CalcLennJonesInteraction_SF( real*, VecR*, real, real, real, real, VecR );
void DistributeVirial( RMat, int, int ); 
void PrintMatrix( RMat );

void ComputeNonBondedInteractions( real *systemPotentialEnergy, /* input/output -- potential energy of the entire system */
				   RMat *systemVirial ) /* input/output -- virial of the entire system */
{
  int i, j, n; 
  real rr, ijPotEn;
  RMat ijVirial, ijVirialSum;  /* used in the accumulation of the total system virial */
  VecR ijForce, jiForce, rij, rChgDip, ijTorqueTot, jiTorqueTot, ijForceTot, ijTorque, jiTorque, torque;
  real potEnLennJonesSum = 0., potEnWatSum = 0., potEnChgSum = 0., potEnChgDipSum = 0., potEnDipSum = 0.;
  real watDipoleSq = Sqr( watDipole );
  real sigWatSq = Sqr( sigWat ); 
  real rrCutLipLip = Sqr( rCutLipLip );
  real rrCutWatWat = Sqr( rCutWatWat );
  real rrCutMax = Sqr( rCutMax );

  MZero( ijVirialSum.u );

  for ( n = 0; n < nebrTabLen; n++ ) {
    i = nebrTab[ 2 * n ];
    j = nebrTab[ 2 * n + 1 ];

    VSub( rij, siteCoords[ i ], siteCoords[ j ] ); 
    VWrapAll( rij );   /* Minimum Image Convention */
    rr = VLenSq( rij ); 
    if ( rr < rrCutMax ) { 
      if ( interactionType[ n ] == WAT_WAT ) { /* water-water interaction (special case) */
	if ( rr < rrCutWatWat ) {  
	  //	  CalcDipDipInteraction(&ijPotEn,&ijForceTot,&ijTorqueTot,&jiTorqueTot,watDipoleSq,siteOrientVec[i],siteOrientVec[j],rij,rr);
	  //	  CalcDipDipInteraction_CS( &ijPotEn, &ijForceTot, &ijTorqueTot, &jiTorqueTot, watDipoleSq, siteOrientVec[i], siteOrientVec[j], rCutWatWat, rij, rr );
	  CalcDipDipInteraction_SF( &ijPotEn, &ijForceTot, &ijTorqueTot, &jiTorqueTot, watDipoleSq, siteOrientVec[i], siteOrientVec[j], rCutWatWat, rij, rr );
	  potEnWatSum += ijPotEn;
	  //	  CalcLennJonesInteraction( &ijPotEn, &ijForce, sigWatSq, epsWat, rr, rij );
	  CalcLennJonesInteraction_SF( &ijPotEn, &ijForce, sigWatSq, epsWat, rr, rrCutWatWat, rij );
	  potEnWatSum += ijPotEn;
	  VVAdd( ijForceTot, ijForce ); 
	  VOuter( ijVirial.u, rij, ijForceTot );
	  if ( (lpp) && ( stepCount % stepLpp == 0 ) ) DistributeVirial( ijVirial, i, j ); 
	  MMAdd( ijVirialSum.u, ijVirial.u );
	  VVAdd( siteForce[ i ], ijForceTot );
	  VVSub( siteForce[ j ], ijForceTot ); /* Newton's Third Law */
	  VVAdd( siteTorque[ i ], ijTorqueTot );
	  VVAdd( siteTorque[ j ], jiTorqueTot );
	} /* end if ( rr < rrCutWatWat ) */	
      } /* end if ( interactionType[ n ] == WAT_WAT ) */
      else { /* evaluate all remaining (non ssd-ssd) interactions */
	if ( rr < rrCutLipLip ) { 
	  CalcLennJonesInteraction_SF( &ijPotEn, &ijForceTot, sigSquared[site[i].type][site[j].type], 
				       eps[site[i].type][site[j].type], rr, rrCutLipLip, rij );
	  potEnLennJonesSum += ijPotEn;
	  switch ( interactionType[ n ] ) { /* possibly evaluate electrostatic interactions */
	  case CHG_DIP: /* iCharge - jDipole */
	    CalcChgDipInteraction_SF( &ijPotEn, &ijForce, &torque, productElectroMoments[site[i].type][site[j].type], 
				      siteOrientVec[j], rCutLipLip, rij, rr );
	    potEnChgDipSum += ijPotEn;
	    VVAdd( ijForceTot, ijForce );
	    VVAdd( siteTorque[ j ], torque );
	    break;
	  case DIP_CHG: /* iDipole - jCharge: the charge-dipole interaction is really evaluated, with the appropriate input/output variables  */
	    VSCopy( rChgDip, -1., rij ); /* input distance vector needs to be rChgDip = rChg - rDip = -rij */
	    CalcChgDipInteraction_SF( &ijPotEn, &jiForce, &torque, productElectroMoments[site[i].type][site[j].type], 
				      siteOrientVec[i], rCutLipLip, rChgDip, rr );
	    potEnChgDipSum += ijPotEn;
	    VVSub( ijForceTot, jiForce ); /* subtract because "jiForce" is the force on the charge site but we need that on the dipole site */ 
	    VVAdd( siteTorque[ i ], torque ); 
	    break;
	  case DIP_DIP: /* dipole-dipole */ 
/* 	    CalcDipDipInteraction_CS( &ijPotEn, &ijForce, &ijTorque, &jiTorque, productElectroMoments[site[i].type][site[j].type], */
/* 				      siteOrientVec[i], siteOrientVec[j], rCutLipLip, rij, rr ); */
	    CalcDipDipInteraction_SF( &ijPotEn, &ijForce, &ijTorque, &jiTorque, productElectroMoments[site[i].type][site[j].type],
				      siteOrientVec[i], siteOrientVec[j], rCutLipLip, rij, rr );
	    potEnDipSum += ijPotEn;
	    VVAdd( ijForceTot, ijForce );
	    VVAdd( siteTorque[ i ], ijTorque );
	    VVAdd( siteTorque[ j ], jiTorque );
	    break;
	  case CHG_CHG: /* charge-charge */
	    CalcChgChgInteraction_SF( &ijPotEn, &ijForce, productElectroMoments[site[i].type][site[j].type], rCutLipLip, rr, rij );
	    potEnChgSum += ijPotEn;
	    VVAdd( ijForceTot, ijForce ); 
	    break;
	  }
	  VOuter( ijVirial.u, rij, ijForceTot );
	  if ( (lpp) && ( stepCount % stepLpp == 0 ) ) DistributeVirial( ijVirial, i, j );
	  MMAdd( ijVirialSum.u, ijVirial.u );
	  VVAdd( siteForce[ i ], ijForceTot );
	  VVSub( siteForce[ j ], ijForceTot );  /* Newton's Third Law... */
	} /* end if ( rr < rrCutLipLip ) */
      } /* end else */
    } /* end if ( rr < rrCutMax ) */
  } /* end for ( n = 0; n < nebrTabLen; n++ ) */
  *systemPotentialEnergy += potEnLennJonesSum + potEnWatSum + potEnChgSum + potEnDipSum + potEnChgDipSum;
  MMAdd( systemVirial->u, ijVirialSum.u );
  /*   printf("###Virial from nonbonded interactions:\n"); */
  /*   PrintMatrix( ijVirial ); */
#ifdef DEBUG
  printf("potEnSum = %g ", potEnLennJonesSum + potEnWatSum + potEnChgSum + potEnDipSum + potEnChgDipSum ); 
  printf("( potEnLJ = %g, potEnChgDipSum = %g ", potEnLennJonesSum, potEnChgDipSum ); 
  printf("potEnChgSum = %g, potEnDipSum = %g, potEnWatSum = %g)\n",  potEnChgSum, potEnDipSum, potEnWatSum  ); 
  fflush(stdout);
#endif
}

